/*! 
 ***************************************************************************
 * \file BluckleyLeverett.cpp
 * Single phase in porous media in 1D.
 *
 * \par Description
 * This code solves the 1D single phase equation for pressure:
 * \f$ \frac{\phi \mu c}{k}\frac{\partial p}{\partial t} = \nabla^2 p \f$.
 * This equation is solved in the interval \f$ x \in [0,100 cm] \f$.
 * The data used are: 
 * The initial condition for the temperature is \f$ P = 0 \f$
 * and the boundary conditions:  \f$ P_L=2 atm \f$ and \f$ P_R=1 atm\f$.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./singlePhase  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Mon Mar 17 08:47:37 GMT 2008 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
//#include "Storage/DiagonalMatrix.hpp" 
#include "Equations/TwoPhaseEquation.hpp"
#include "./TP/BLIP1.hpp"
#include "./TP/BLES1.hpp"
#include "Solvers/TDMA.hpp"
#include "Solvers/explicit.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"

#include "DDM/CartComm.hpp"
#include "DDM/SubDomain.hpp"

using namespace Tuna;
  
// Scalar field in two dimensions
typedef TunaArray<double,2>::huge ScalarField2D;

timer crono;

//-------------------------------------------------------
// Parameters to make the partition
double lx_g, ly_g;
int num_nodes_x_g, num_vols_x_g, num_subdom_x;
int num_nodes_y_g, num_vols_y_g, num_subdom_y;
int nc_ovlp_l, nc_ovlp_r, nc_ovlp_d, nc_ovlp_u;
int Isub, Jsub, rank, size;

// Local variables for each subdomain
int num_nodes_x, num_nodes_y;
double length_x, length_y;
double dt, tolerance = 1e-6, Tmax;
int tdma_iter = 20, steps, frequency;
double porosity, permeability, mu_w, mu_o, Srw, Sro;
double injection, pres_out, Sw_in;
int idt_p;

bool read_data_and_Bcast(CartComm<2>&);
//-------------------------------------------------------
int main( int argc, char **argv)
{


    MPI::Init(argc, argv);
    rank = MPI::COMM_WORLD.Get_rank();       
    size = MPI::COMM_WORLD.Get_size(); 
    
// ----- Construction of a Cartesian topolog
    CartComm<2> cart(argc, argv, rank);   
    Isub = cart.get_I();
    Jsub = cart.get_J();
    num_subdom_x = cart.getNumProc_I(); 
    num_subdom_y = cart.getNumProc_J();
    //    cart.print();    

// Reading and broadcasting data from the input file 
    read_data_and_Bcast(cart);   

// Global Mesh
    StructuredMesh<Uniform<double, 2> > mesh(lx_g, num_nodes_x_g,
					     ly_g, num_nodes_y_g);
    double dx = mesh.getDelta(X);
    double dy = mesh.getDelta(Y);

    //    mesh.print();

    /***/
// Subdomain
    SubDomain<double, 2> subdom(cart);
    
    double ovlp_l = subdom.createOverlap(LEFT,  nc_ovlp_l, dx);
    double ovlp_r = subdom.createOverlap(RIGHT, nc_ovlp_r, dx);
    double ovlp_d = subdom.createOverlap(DOWN,  nc_ovlp_d, dy);
    double ovlp_u = subdom.createOverlap(UP,    nc_ovlp_u, dy);
  
// Local mesh for subdomains of equal length
    length_x = lx_g / num_subdom_x; // Non-overlapping
    length_y = ly_g / num_subdom_y; // Non-overlapping
    
    // Restriction: (num_nodes_x_g - 1) % num_subdomx = 0 
    num_nodes_x = (num_nodes_x_g - 1) / num_subdom_x + 1; 
    // Restriction: (num_nodes_y_g - 1) % num_subdomy = 0 
    num_nodes_y = (num_nodes_y_g - 1) / num_subdom_y + 1; 

// Adding the overlapping
    length_x += ovlp_l + ovlp_r; 
    length_y += ovlp_d + ovlp_u; 
    num_nodes_x += nc_ovlp_l + nc_ovlp_r;
    num_nodes_y += nc_ovlp_d + nc_ovlp_u;

// Resizing the mesh to the local values for the subdomain
    mesh.resize(length_x,  num_nodes_x, length_y, num_nodes_y);

    mesh.print();
// Define the info to be exchanged between subdomains
    subdom.setupInfoExchange(mesh);
    subdom.print();

    double Gamma = permeability;

    ScalarField2D p ( mesh.getExtentVolumes() );
    ScalarField2D Sw ( mesh.getExtentVolumes() ); 
    ScalarField2D p_n( mesh.getExtentNodes()   );
    ScalarField2D Sw_n( mesh.getExtentNodes()   );

    // Initial conditions
    Sw = 0;
    p = pres_out; 
 
    Range all = Range::all();
    // Boundary condition
    Sw(0, all) = Sw_in;  

    /** */
    NumUtils::interpolateToNodes(p_n, p);
    NumUtils::interpolateToNodes(Sw_n, Sw);

    if (rank == 0) {
      InOut::writeToFile_DX(Sw_n, 0, "./DataBL/00/satu.", dx, dy);
      InOut::writeToFile_DX(p_n, 0, "./DataBL/00/pres.", dx, dy);    
    }    
    if (rank == 1) {
      InOut::writeToFile_DX(p_n, 0, "./DataBL/01/pres.", dx, dy);    
      InOut::writeToFile_DX(Sw_n, 0, "./DataBL/01/satu.", dx, dy);
    }
    if (rank == 2) {
      InOut::writeToFile_DX(Sw_n, 0, "./DataBL/02/satu.", dx, dy);
      InOut::writeToFile_DX(p_n, 0, "./DataBL/02/pres.", dx, dy);    
    }    
    if (rank == 3) {
      InOut::writeToFile_DX(p_n, 0, "./DataBL/03/pres.", dx, dy);    
      InOut::writeToFile_DX(Sw_n, 0, "./DataBL/03/satu.", dx, dy);
    }
    
    //    SparseMatrix< Diagonal< double, 2> > A(num_nodes_x, num_nodes_y); 
    DiagonalMatrix< double, 2> A(num_nodes_x, num_nodes_y); 
    ScalarField2D              b(num_nodes_x, num_nodes_y);

    TwoPhaseEquation< BLIP1<double, 2> > pressure(p, A, b, mesh.getDeltas());
    pressure.setDeltaTime(dt); 
    pressure.setPermeability(permeability);
    pressure.setPorosity(porosity);
    pressure.setSrw(Srw);
    pressure.setSro(Sro);
    pressure.setViscosity_w(mu_w);
    pressure.setViscosity_o(mu_o);
    pressure.setNeumann  (LEFT_WALL, -injection * mu_w / permeability);
    pressure.setDirichlet(RIGHT_WALL);
    pressure.setNeumann  (TOP_WALL);
    pressure.setNeumann  (BOTTOM_WALL);
    pressure.setSaturation(Sw);
    pressure.print();

    //
    // In these case, the boundary conditions are calculated in
    // the numerical scheme (the adaptor BLES1).
    //
    TwoPhaseEquation< BLES1<double, 2> > saturation(Sw, A, b, mesh.getDeltas());
    saturation.setDeltaTime(dt); 
    saturation.setPermeability(permeability);
    saturation.setPorosity(porosity);
    saturation.setSrw(Srw);
    saturation.setSro(Sro);
    saturation.setViscosity_w(mu_w);  
    saturation.setViscosity_o(mu_o); 
    saturation.applyBounds(1, Sw.ubound(firstDim)-1);// just to move the indexes
    saturation.setNeumann(TOP_WALL);
    saturation.setNeumann(BOTTOM_WALL);
    saturation.setPressure(p);
    saturation.print();

    ofstream error_fp, residual_fp;
    
    if (rank == 0){
      error_fp.open("./DataBL/00/error"); 
      residual_fp.open ("./DataBL/00/residual");
    }
    if (rank == 1) {
      error_fp.open("./DataBL/01/error"); 
      residual_fp.open("./DataBL/01/residual");
    }
 
    int iteration = 1, idprint = 1;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    //    int pausa;

    int bi = saturation.get_bi(), bj = saturation.get_bj();
    int ei = saturation.get_ei(), ej = saturation.get_ej();

    frequency *= (86400 / dt);

    double total_time = 0, pres_time = 0, sat_time = 0;
    while (t <= Tmax) {

      if ( !(iteration % idt_p) ) {
	crono.tic();
	pressure.calcCoefficients();
	Solver::TDMA2DX(pressure, tolerance, tdma_iter, 1.0);
	pres_time += crono.toc();
	
	residual = pressure.calcResidual();
	error = pressure.calcErrorL1();
	
	error_fp << iteration << "\t" << error << std::endl;
	residual_fp << iteration << "\t" << residual << std::endl;	
	
	pressure.update();
      }  

      crono.tic();
      saturation.calcCoefficients();
      Solver::solExplicit2D(saturation);
      sat_time += crono.toc();

      error = saturation.calcErrorL1();
      saturation.update();



      /* */
      // just to avoid go out of boundary conditions.
      for(int j = 0; j <= (ej+1); ++j)
	for(int i = 0; i <= (ei+0); ++i) {
	  if ( Sw(i,j) > Sw_in) Sw(i,j) = Sw_in;
	  if ( Sw(i,j) <= 0 ) {
	    Sw(i+1,j) = 0.0;
	  }
	}
      saturation.updatePhi(Sw);
      /* */
     
      if ( !(iteration % frequency) ) {

	NumUtils::interpolateToNodes(p_n, p);
	NumUtils::interpolateToNodes(Sw_n, Sw);

	if (rank == 0) {
	  InOut::writeToFile_DX(p_n, idprint, "./DataBL/00/pres.",dx, dy);	
	  InOut::writeToFile_DX(Sw_n, idprint, "./DataBL/00/satu.",dx, dy);
	}
	if (rank == 1) {
	  InOut::writeToFile_DX(p_n, idprint, "./DataBL/01/pres.",dx, dy);	
	  InOut::writeToFile_DX(Sw_n, idprint, "./DataBL/01/satu.",dx, dy);
	}
	if (rank == 2) {
	  InOut::writeToFile_DX(p_n, idprint, "./DataBL/02/pres.",dx, dy);	
	  InOut::writeToFile_DX(Sw_n, idprint, "./DataBL/02/satu.",dx, dy);
	}
	if (rank == 3) {
	  InOut::writeToFile_DX(p_n, idprint, "./DataBL/03/pres.",dx, dy);	
	  InOut::writeToFile_DX(Sw_n, idprint, "./DataBL/03/satu.",dx, dy);
	}

	std::cout << "\n\n ---- ID print => " << idprint
		  << "\t - Iteration => " << iteration
		  << "\t Step = " << t
		  << "\n ---- Res = " << residual
		  << "\t Err = " << error;
	//std::cout.flush();
	++idprint;
      }
    

      /* */
      subdom.infoExchange1(p);
      pressure.updatePhi(p);

      subdom.infoExchange1(Sw);
      saturation.updatePhi(Sw);
      /* */
      // just to avoid go out of boundary conditions.
      for(int j = 0; j <= (ej+1); ++j)
	for(int i = 0; i <= (ei+0); ++i) {
	  if ( Sw(i,j) > Sw_in) Sw(i,j) = Sw_in;
	  if ( Sw(i,j) <= 0 ) {
	    Sw(i+1,j) = 0.0;
	  }
	}
      saturation.updatePhi(Sw);
      /* */

      ++iteration;
      t += dt;
    }

    cout << "\n Pres Etime = " << pres_time << " segs \n";
    cout << "\n Sat Etime = " << sat_time << " segs \n";

/* **/

    MPI::Finalize();
    return 0;
}





bool read_data_and_Bcast(CartComm<2>& cart) 
{
  int    data_Int[9];
  double data_Flt[12];


  if (rank == 0) {
    std::cout << "\n\n"
	      << "          MASTER PROCESS rank = " << rank << "\n"
	      << " +----------------------------------------+\n"
	      << " |       TUNA::FVM FOR PDE SOLVING        |\n"
	      << " +----------------------------------------+\n"
	      << " | No subdomains in x-axis " << num_subdom_x << "\n"
	      << " | No subdomains in y-axis " << num_subdom_y << "\n"
	      << " +----------------------------------------+\n"
	      << "\n";
	
    // ----- Reading data from "input" file	
    std::ifstream input_cin ("inputBL2DPara");
    input_cin >> data_Flt[0]  // length in x-axis
	      >> data_Flt[1]  // length in y-axis
	      >> data_Int[0]  // Num of nodes in x-axis for each subdomain
	      >> data_Int[1]  // Num of nodes in y-axis for each subdomain
	      >> data_Int[2]  // Overlapping left
	      >> data_Int[3]  // Overlapping right
	      >> data_Int[4]  // Overlapping down
	      >> data_Int[5]  // Overlapping up
	      >> data_Flt[2]  // Porosity
	      >> data_Flt[3]  // Permeability
	      >> data_Flt[4]  // mu_w
	      >> data_Flt[5]  // mu_o
	      >> data_Flt[6]  // Srw
	      >> data_Flt[7]  // Sro
	      >> data_Flt[8]  // Injection
	      >> data_Flt[9]  // pres_out
	      >> data_Flt[10] // Sw_in
	      >> data_Flt[11] // dt Time step
	      >> data_Int[6]  // Tmax Time steps
	      >> data_Int[7]  // idt_p
	      >> data_Int[8];  // Print frequency
     
    input_cin.close();
  }
// ----- Broadcast the info to all processors
  cart.comm.Bcast(data_Int, 9, MPI::INT, 0);    
  cart.comm.Bcast(data_Flt, 12, MPI::DOUBLE, 0); 

// ----- Using global variables.
    lx_g    = data_Flt[0];
    ly_g    = data_Flt[1];

    num_nodes_x_g = data_Int[0];
    num_nodes_y_g = data_Int[1];
    nc_ovlp_l     = data_Int[2]; 
    nc_ovlp_r     = data_Int[3];
    nc_ovlp_d     = data_Int[4]; 
    nc_ovlp_u     = data_Int[5];

    porosity     = data_Flt[2];
    permeability = data_Flt[3];
    mu_w         = data_Flt[4];
    mu_o         = data_Flt[5];
    Srw          = data_Flt[6];
    Sro          = data_Flt[7];
    injection    = data_Flt[8];
    pres_out     = data_Flt[9];
    Sw_in        = data_Flt[10];
    dt           = data_Flt[11];

    Tmax         = data_Int[6];
    idt_p        = data_Int[7];
    frequency    = data_Int[8];


    return 0;
}


